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ABSTRACT 

Aims. Cold steady-state disk wind theory from near Keplerian accretion disks requires a large scale magnetic field at near equiparti- 
tion strength. However the minimum magnetization has never been tested with time dependent simulations. We investigate the time 
evolution of a Shakura-Sunyaev accretion disk threaded by a weak vertical magnetic field. The strength of the field is such that the 
disk magnetization falls off rapidly with radius. 

Methods. Four 2.5D numerical simulations of viscous resistive accretion disk are performed using the magnetohydrodynamic code 
(*y| PLUTO. In these simulations, a mean field approach is used and turbulence is assumed to give rise to anomalous transport coefficients 

, (alpha prescription). 

Results. The large scale magnetic field introduces only a small perturbation to the disk structure, with accretion driven by the dom- 
inant viscous torque. However, a super fast magnetosonic jet is observed to be launched from the innermost regions and remains 
stationary over more than 953 Keplerian orbits. This is the longest accretion-ejection simulation ever carried out. The self-confined 
jet is launched from a finite radial zone in the disk which remains constant over time. Ejection is made possible because the mag- 
netization reaches unity at the disk surface, due to the steep density decrease. However, no ejection is reported when the midplane 
magnetization becomes too small. The asymptotic jet velocity remains nevertheless too low to explain observed jets. This is because 
' of the negligible power carried away by the jet. 

Conclusions. Astrophysical disks with superheated surface layers could drive analogous outflows even if their midplane magnetiza- 
tion is low. Sufficient angular momentum would be extracted by the turbulent viscosity to allow the accretion process to continue. 
The magnetized outflows would be no more than byproducts, rather than a fundamental driver of accretion. However, if the midplane 
magnetization increases towards the center, a natural transition to an inner jet dominated disk could be achieved. 

Key words, accretion, accretion disks - Magnetohydrodynamic s (MHD) - stars: formation - ISM: jets and outflows - galaxies: 
nuclei - galaxies: jets 



Q ■ 1. Introduction iBlandford & Pavnd (Il982l) establishes a relationship between 

' . ,. , , . . the mass loading and the magnetic lever arm of magnetocen- 

^ ■ Accretion disks are commonly found in young stars active trif n driven outflows But the etic field st th was 

galactic nuclei, cataclysmic variables and microquasars. In order kft unconstrainedi so in principal magn etization at the disk 

to allow material to accrete onto a central object, it is necessary surface couM driye a low . enthal outflow xhe reason lies in 

^ ; to lose some angular momentum in an efficient way. This is pos- ^ ^ ^ an ideal MHD jgt modd assumefj ^ mafjs losfj 

. sible in a disk in one of two ways, eith er by radial outward trans - and doef , QQt ufc k afj function of ^ digk eters xhis 



port in a ^disk by turbulent tra nsport QMssUl Sunyaev 19731 was precise i y the g oal of semi -analyti cal studies done by e.g. 
Lvnd en-Bell& Pringlel | l974 or spiral waves jTagger & Pellat| | pelletier ( m ^ and gf£y (J^ where ^ digk 



r = — : — , " * 1 , _ ; — rr~—. r rerrerra 6i renetier ( ivvd) ana rerreira (ivy/), wnere tne aisK 

999)^r_j^aLtes|^rt upwards out of the disk in a jet h fe consist ently computed: these authors showed 

(Blandford& Payne 1982). - ! 



that steady-state cold jets can be produced only with a vertical 



Two extreme possible disk structures can then be identified, fi ij i , • *f „ • ,., • t f f . 

r ue\a close to equipartition. A few numerical experiments tested 

corresponding to each of these two processes of angular momen- ... • ■ . * 

F 6 F 5 the accretion-ejection connection in a consistent way: axisym- 



tum removal. 



„ , _ . . , ... . metric magneto-hydrodynamic (MHD) simulations of resistive 

The Jet Emitting Disk (hereafter JED) is threaded by a large accretion disks rting ^ produ ct i n of se lf-confined, quasi- 

scale magnetic field of bipolar topology driving a jet (defined stead super . fast jets (Casse & KeTOens 2002i 2004; Zanni et al . 

here as super-fast magnetosonic flow). The dominant torque 2007; Xzeferacos et al 2009) confirmed most of the resu i ts ob _ 

m the JED is magnetic, due to the large braking lever arm of tained ^ semi . anal tical models . The were however done 
the jet, defined by a length scale equivalent to the Altven ra- 
dius (Pelletier & Pudritzl ll992h . The pioneering jet model by 



Send offprint requests to: G. C. Murphy 



2 



Gareth C. Murphy et al.: Ejection from weakly magnetized disk 



with a large disk magnetization^ in the range fi ~ 0.2 - 1 . The 
inner regions of the disk from whence jets are observed to be 
emitted are expected to be JED-like. In the specific case of out- 
flows from young stars, extrapolation of slitless images of Class 
II jets have constrained the launching region to be confined to a 
zone of radial extent ~ a few AU close to the centre of the disk 
(iHartigan et al.ll2Q04t ICabritll2Q07h . 

The outer regions are expected to behave more like 
the well studied standard accretion disk , hereafter SAD 
( Shakura & Sunvaevlll97l iFrank et al.ll2002h . Here the charac- 
teristic lengthscale over which the viscous torque is exerted is 
of the order of ~ a w h, where a v measures the level of the tur- 
bulence. Such a turbulence is assumed to arise from the devel- 
opment of magnetic instabilities that are triggered in the dis k 
whenever a magnetic field is present (Bal bus & Hawlevlll991l) . 
This field must however be below equipartition strength (namely 
B 2 //V() <K P) to avoid the stabilizing effect of the magnetic ten- 
sion. Therefore, the high magnetization required by a JED im- 
pedes the development of disk turbulence which is however re- 
quired to support a steady launching: that would leave only a 
very tiny parameter space for stationary ejection to take place. 
The SAD- JED structure h a s been put forward in sever a l papers 



Ferreira et al] d2006l) ; ICombet & Ferreiral (120081) : Eerreira 



(120081) . 

The study of low magnetization accretion regimes has been 
attempted making use of fully 3D global simulations of accre- 
tion disks, threaded by a weak large scale mag netic field. MRI 
sets i n and accretion is quickly established (lHawlev & Balbusl 
|2002|). The remarkable resu lt is that outflows are also produced 
(Igumens hchev et al]|2003[). especially when the imposed field 
is of bipolar topology ( Beckwith et al.l l2009l) . However, many 
questions remain open: What controls the mass loss in these sim- 
ulations? Will the outflowing plasma become a self-confined jet? 
Is grid resolution enough to properly describe the turbulent cas- 
cade? The fact is that it is still impossible to properly follow tur- 
bulence while solving for the long term evolution of large scale 
systems. 

As a consequence, the question of super-fast magnetosonic 
jet formation from weakly magnetized disks is still open. In 
this paper we address this issue using 2.5D numerical MHD 
simulations based on a mean field approximation. We explore 
the accretion-ejection processes from a quasi-standard accre- 
tion disk where the magnetization is very low (smaller than 
10~ 3 ). Since the magnetic field is low, we assume that turbu- 
lence triggered by the MRI is indeed present but that it pro- 
vides mainly anomalous transport coefficients: a viscosity v v 
and a magnetic diffusivity v m . On the other hand, we do not 
expect to observe any MRI feature (such as channel flows for 
instance) in our simulation because of the presence of explicit 
viscosity and magnetic diffusivity effects. While measurements 
of the turbulent viscosity in MRI induced turbulence have been 
extensively reported in the literature, it is only very recently that 
such a work has been done for the turbulent magnetic diffusivity 
dLesur & Longarettil2009l:lGuan & Gammidl2009l) . In particular 
Lesur & Longaretti (2009) showed that the turbulent magnetic 
diffusion scales like a resistivity tensor with dominant diagonal 
terms. Also, as a first approximation, an isotropic value can be 
safely used. Finally, the effective Prandtl number V m = v v /v m , 
given by the ratio of turbulent viscosity and diffusivity, has been 
found to be of order unity. The mean field approximation has 



1 The magnetization is related to the usual plasma beta by \i = 2//? in 
gas pressure supported disks. It is however a more general concept as s 



been successfully employe d in a number of semi-analytical (e.g . 
i Ferreira & Pelletieri 119951: ID 119951: ICasse & Ferreiral |2000a; 
lOgilvie & Livid 1200 U |Rothstein & Lovelace 2008 ) and nu- 
merical applications ("e.g. JCasse & Keppens 2002; Kukeretal 
2003; von RekowskTet alj2003: Melia ni et alj2006HZanni et al 



2007; iRomanova et al.l l2009) related to the study of magnetized 
accretion-ejection flows. Beside having a precise control of the 
diffusive and transport phenomena, the numerical experiments 
based on this approach provide laminar flow solutions which can 
be compared to semi-analytical models. 

In section 2, we describe the numerical method used, the 
boundary and initial conditions. Section 3 is devoted to the de- 
scription and discussion of the results obtained. Surprisingly, 
super-fast jets are indeed obtained from a finite disk region and 
remain stable for a time span never previously achieved in the lit- 
erature. Section 4 summarizes our findings and, in a companion 
paper (Murphy et al., in prep), we will examine the long stand- 
ing issue of the magnetic field redistribution within the disk on 
long (accretion) time scales. 

2. Numerical method 

The full visco-resistive MHD equations in axial symmetry are 
evolved i n time using the pu blicly available numerical code 
PLUTO (Mi gnone et alj|2007h . The solved equations are: the 
continuity equation 



dp 
dt 



+ V • (pu) = : 



(1) 



the conservation of momentum equation 



pu®u + CP*) I + B®B + T 
the induction equation 



+ p VO G = ; (2) 



dB 

— + Vx(Bxh + v m J) = ; 
at 



the conservation of energy equation 

(E + P*)u-{u- B)B + v m J xB-u T 



9E TT 

dt 



= S 



(3) 



(4) 



where S = -puV<S>c + L c , and L c is the local cooling term (see 
below). The total energy density is defined as 



E = \p\u\ 2 + JL- + \\B\ 2 , 
2 y— 1 2 

and the total pressure (thermal and magnetic) is 

P* =P+-\B\ 2 . 
2 



(5) 



(6) 



The equations are written and solved in dimensionless form, 
thus without juo coefficients. The equation of state is the ideal 
gas equation. Here, p is the mass density, u the velocity, P the 
gas pressure, B the magnetic field, <&g = -GM/ Vr 2 + z 2 is the 
gravitational potential of the central mass, J = V x B is the 
current density, v m the magnetic diffusivity and y - 5/3 is the 

ratio of specific heats. The viscous stress tensor T is defined as 



it is defined with the total pressure P gas + P r 



ad- 



(Vm) + (Vm) 7 



3 v ' 



(7) 
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where rj v is the dynamic viscosity. See Appendix lAl for the ex- 
pression of the tensor components. As is customary, the kine- 
matic viscosity is defined as v v = J] v /p. 

As stressed above, we follow a mean field approach where 
the turbulence is crudely modeled by mere transport coefficients: 
a viscosity v v and a magnetic diffusivity v m . Co nsistently with 
this approximation, a Shakura & Sunvaevl d 19731) alpha prescrip- 
tion is then employed. This assumes that the viscosity is propor- 
tional to the heightscale of the disk, h, and some characteristic 
velocity, in this case the sound speed, c s , namely 



3 v v 
2 c s h 



(8) 



We assume that the disk is not flat, but will have initially a con- 
stant aspect ratio e = h/r = c s /Vk = 0.1. As an initial condition 
for the alpha accretion disk, we take th e perturbative solution of 
the steady-state disk equations found in lZanni & Ferreiral (T2009) 
and the references therein. The disk is in hydrostatic equilibrium 
and accretion is driven by the viscous stress tensor alone. This 
particular solution provides reasonable vertical and radial pro- 
files of all quantities that are suitable for a SAD (see Appendix 
|A]for more details). A not so well known bias of the alpha pre- 
scription in 2D flows is that, below a critical value found to be 
ff cr it ~ 0.685, there is a backflow on the disk midplane (Urpin 
1984) . This is certainly unphysical and arises only from the func- 
tional form of the stress tensor used to mimic turbulence. In or- 
der to circumvent this bias, we u sed a v — O.S0. 

Consistently with the recent iLesur & L ongarettil d2009l) re- 
sults, we assume that the effective magnetic Prandtl number 
f m = v v /v m is of order unity: for simplicity we set V m = 
2/3 in all simulations. Again, we stress that this is a strong 
simplification of highly complex phenomena but also an un- 
avoidable price to pay if one seeks for long term evolution of 
global systems, such as accretion disks and their related jets. 
With a constant !P m , the viscosity and resistivity will follow 
the same radial and vertical profiles. They decrease smoothly 
with height until they become negligible, allowing a transition 
to a magnetized "corona" in ideal MHD regime. Since MRI in- 
duced turbulence is qu enched when the magnet ic field becomes 
close to equipartition (Balbus & Hawlevl [l991l) . there will be a 
height where the accretion flow cannot be turbulent anymore 
(Roth stein & Lovelacel l2008). For simplicity, we assume that it 
corr esponds to the disk surface (see App endix IA1. W e thus f ol- 
low ICasse & Keppensl d2002l |2004l) and IZanniet all d2007l) in 
neglecting the turbulent viscosity and turbulent resistivity in the 
highly magnetised corona. 

In a real accretion disk, the local heating due to turbulence 
(here crudely modeled by alpha prescriptions for resistivity and 
viscosity) would be balanced by both turbulent transport and ra- 
diative cooling. While the former cooling term needs full 3D 
calculations, the latter can be done in 2D but requires radia- 
tive transfer. Both effects are far beyond the scope of the present 
work. Hence, by including a "cooling" function L c such that 

L c = v m J 2 + ^- [T 2 rr + Tl + T 2 ^ + 2(T 2 Z + T% + T%j\ , (9) 

we can exactly balance both resistive and viscous heating terms. 
Our disk evolution is therefore adiabatic despite the presence of 
transport coefficients within the disk. This is certainly a caveat, 



2 This value might be seen too large when compared to the small 
mean field used in the disk. However, note that the main effect of a 
large a v is to reduce the accretion time scale, while still maintaining it 
well below the Keplerian one. 



shared by most today MHD simulations, and deserves further in- 
vestigation. On the other hand, it allows to avoid in a simple way 
the otherwise unavoidable (and unphysical) vertical expansion 
of the disk when heating is present without any kind of cooling. 
A static atmosphere in pressure equilibrium is set above the disk 
and a large scale magnetic field is superimposed in the whole 
domain. 

To set the initial magnetic field, we use the magnetic flux 
function such that B = xe^/r. We take the particular form 



m{r,z)=AB^\- 



1/4 



„7/4 



(m 2 + z 2 /r 2 ) 



(10) 



where Bo = Vw^o^do, Pdo is the thermal pressure of the 
disk and fi is the disk magnetization at the disk midplane of the 
inner radius ro. The parameter m describes the initial bending of 
the magnetic field lines, and is set to 0.935 in all simulations. 
This leads to a magnetic field such that the disk magnetization 
li varies initially as r _1 , starting from //(ro) = 2 x 10~ 3 (one 
simulation is done with /i(ro) = 2 x 10~ 4 ). Such a small value 
for the large scale field is chosen to ensure that it is only a tiny 
perturbation to the initial SAD structure on the midplane. 

The MHD system of equations has been solved numerically 
exploiting the MHD module provided with PLUTO. The code 
has been configured to perform second-order piecewise linear 
reconstruction of primitive variables, with a Van Leer limiter for 
the density and magnetic field components and a minmod lim- 
iter for the thermal pressure and velocity components. To com- 
pute th e intercell fluxes, a HLL Riemann solver has been em- 
ployed dHartenetal.1119831) . while second order in time has been 
achieved using a Runge-Kutta scheme. The solenoidal condition, 
V • B = 0, is preserved using the constrained transport method 
dEvans & Hawlevl [l988). The viscous and resistive terms have 
been treated explicitly, using a second-order finite difference ap- 
proximation for the dissipative fluxes and checking the diffusive 
timestep. 

A uniform resolution grid of 512 by 1536 cells is used. This 
describes a domain of 40ro by 120ro, where ro is the inner radius 
of the disk. An outer, stretched grid is extended for a further 512 
cells in the radial direction and 1536 in the vertical direction, 
thus describing in total a region 280ro by 840ro. To examine the 
effects of a higher resolution (see Section [3~5l l. the same grid is 
used again but this time only describing a region of 5ro by 15ro, 
with a logarithmically stretched grid outside this region 35ro by 
105ro. A disadvantage of higher resolution is that the number of 
timesteps to reach accretion timescales becomes prohibitive. 

The boundary conditions are axial symmetry on the rota- 
tion axis and equatorial symmetry for the disk midplane. The 
upper r and z boundaries border on a logarithmically stretched 
grid which ensures that the magnetized outflow never reaches 
the boundaries. The ghost cells at the upper r and z boundaries 
are set to equal the values inside the domain (the numerical ap- 
proximation to an "outflow" boundary condition). The gravita- 
tional potential has a singularity at the origin, so a rectangular 
portion of the simulation close to the origin is excluded, as in 
ICasse & Keppensl J2002 ) . The right boundary of the rectangular 
region is a sink and the upper boundary injects a small amount 
of material into the grid at the escape velocity (with density 1.1 
times the local initial density). This keeps the axis sufficiently 
dense to ensure that unphysically low densities are not produced 
on the axis by the Lorentz force. For the poloidal magnetic field 
the boundary condition can be expressed in terms of the toroidal 
electric field, E$, Assuming flux is not advected into the central 
object, we impose = 0. 
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Throughout the paper the distances are expressed in units of 
ro, which is the inner disk radius in the simulation. Velocities 
have been normalized on the Keplerian speed at ro, Vk,o = 
VGM/ro- Densities are expressed in units of pdo, the initial 
disk density at its inner radius. The times are in units of the 
Keplerian orbital period tko = 277To/Vk,o- Pressures are given 
in units of pdoV^ while the magnetic field is expressed in units 

of ^oPdoV^ . 

For ease of reproducibility, the C subroutines defining initial 
conditions and boundary conditions are available from the au- 
thors on request. The numerical code PLUTO is publicly avail- 
able from the URL http://plutocode.to.astro.it 

3. Ejection from weakly magnetized disks 

3. 1 . Global description 

When the simulations starts in the first visible phenomenon 
is the triggering of the familiar vertical torsional Alfven wave 
dMouschovias & Paleologou 1980; lQuved & Pudritzll 1997b . It is 
due to the differential rotation between the Keplerian disk and 
the initially non rotating atmosphere. But after a few inner disk 
rotations, a proper MHD outflow is launched from the disk, de- 
veloping a bow shock and compressing the ambient material and 
the preceding torsional flow. Figure[T]shows a plot of the jet den- 
sity in the poloidal plane together with the fast and Alfven sur- 
faces and magnetic field lines. A superfast jet is launched within 
a relatively narrow region at the disk surface up to r — 5. Matter 
launched from this region crosses the slow and Alfven surfaces 
close to the disk surface and is accelerated up to the fast magne- 
tosonic surface. 

Along the z direction, the numerical simulation can be char- 
acterised as divided into two main zones, a resistive zone, where 
resistive effects are important (the disk), and an ideal MHD 
zone, where ideal MHD is strictly enforced (the jet and atmo- 
sphere). The disk surface could be defined as the altitude where 
all transport coefficients vanish. We choose rather to define the 
disk surface as the altitude where the radial velocity compo- 
nent vanishes, marking therefore a clear transition between un- 
derlying accreting layers (u r < 0) and outflowing upper layers 
(u r > 0). Figure [2] shows these various surfaces at the final time 
t — 953tko. Note however that they do not evolve much over 
time as it can be seen in Fig.fT] 

This result is unexpected since the normal requirement for 
a steady cold MH D disk wind is a near e quipartition midplane 
magnetization dFerreira & Pelletierf[l995l) . A second surprising 
feature is that the outflow is launched from a clearly defined re- 
gion centrally located in the disk that does not evolve from the 
entire disk region, even after 953 disk rotations. This is in con- 
trast to the previous results of IZanni et al.l d2007l) where, as the 
simulation evolves in time, the outflowing region moves outward 
on the Keplerian timescale. In fact, the global accretion-ejection 
configuration exhibits three distinct zones in the the radial direc- 
tion. Zone I corresponds to the innermost radii where anchored 
field lines give rise to a super-fast jet, namely from r — 1 to 
r = 5. Then an intermediate zone gives birth to a sub-fast but 
still super- Alfvenic outflow. Zone II goes from r = 5 to r = 13. 
Such a zone is expected to be unsteady as any FM wave can 
travel upstream. Finally, the last zone goes from r = 13 up to 
the outermost radius and corresponds to negligible outflowing 
material that remains always sub-Alfvenic. 

This is the longest accretion disk simulation ever done so 
far (953 inner periods) where the jet remains steady. Figure [3] 





It 

disk 


1 








Fig. 2. Poloidal cross section showing various zones at a time 
t=953TKo: the critical surfaces of the MHD outflow (slow mag- 
netosonic Ssm, Alfven Sa, fast magnetosonic Sfm) an( l the disk 
surface. Field lines anchored at r = 1,5.1, 13 are also shown, 
delimiting the three zones (see text) 



shows the accretion (measured on one half disk thickness) and 
ejection (in one jet) rates as well as their ratio plotted over time. 
M w is computed at the disk surface, defined by the height where 
the poloidal velocity reaches zero, and from a radius r = 1 .4 to 
r = 5.0. As can be seen from the figure, the ejection to accretion 
mass loss is approximately 7%. While smaller than that obtained 
for a disk with a larger magnetization, this mass loss is by no 
means negligible. We also inject a small amount of mass M; n j at 
the surface of the internal boundary which is of the order of 1 % 
of M w and thus negligible when compared to M w . 

The accretion power is computed as the difference between 
the flux of mostly mechanical energy E — — + + <Dq en- 
tering the disk at its outer edge and leaving it at its inner edge, 
namely 



I pEUp ■ dS - I 

«7out *Jin 



pEu p 



dS 



(ID 



where the integration is performed on a vertical section of the 
disk. The jet power is calculated as the sum of all energy fluxes 
(mechanical and Poynting) leaving the disk, 



^jet — ^mechjet + ^MHDjet 



where 

-^mechjet 
^MHDjet 



= J pEu p ■ dS 



ExB-dS 



(12) 



(13) 



(14) 



Here, the integration has been made in a control volume defined 
by the inner radius r; n = 1.4 and an outer radius r out = 5.0. 
The theoretical global energy budget should then be 



^acc "t" ^visc — ^jet ^rad 



(15) 



where P m d is the power released into heat by both viscous and 
Joule terms (and is eventually radiated at the disk surface) and 

^visc = J in (u -T) ■ dS is the influx of energy at the inner radius due 
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Fig. 1. Log of mass density at times (a) t = 
0t K0 , (b) t = 31t ko , (c) t = 160t K o, (d) / = 
953tko- The fast surface and the Alfven surface 
are overplotted with dotted and dashed lines 
respectively. The super fast-magnetosonic out- 
flow (the jet) is launched only from a small in- 
ner region, located between r = 1 and r = 5. 
The extension of this zone remains constant 
over time. 



to viscosity (the flux at the outer radius is negligible). In the stan- 
dard accretion disk theory, such a flux of energy is set to exactly 
zero through the "zero torque condition" at the inner boundary. 
This was not implemented in our simulation so that the actual 
power that is available (liberated) within our simulation box is 
^lib = ^acc + ^visc- Also, in practice, our simulation does not in- 
clude radiation but the Joule and viscous heating terms are bal- 
anced exactly by the cooling term L c . Table[TJshows the different 
calculated powers. f v ; sc represents only 20% of all the liberated 
power (Pyisc/^acc = 0.25): while not strictly negligible, it is only 
a small fraction. In the following, we will thus compare the jet 
power only with the accretion power. The total jet power (MHD 
+ mechanical) represents only 15% of the the accretion power. In 
terms of released power, the disk behaves therefore exactly like 
a standard accretion disk, the jet being a mere epiphenomenon. 
The obvious reason for that is the very weak magnetic field, un- 
able to extract a significant fraction of the available power (see 
below). Note also that the jet power is dominated by the MHD 
Poynting flux. 

The jet becomes super-SM very soon, almost at the disk sur- 
face and reaches the Alfven speed at an altitude za significantly 
smaller than the corresponding Alfven radius r& (Fig. [2]). This 
is again in str ong contrast w ith self-similar cold jet solutions 
where za ~ '"a (Ferreira 1997). The flow then reaches its asymp- 
totic velocity soon after the fast magnetosonic surface, which 
is a maximum of about ~ 1.2 times the Keplerian value at the 
disk midplane (see Fig. |4|. Thus the type of jet produced here 
cannot be responsible for very high velocity Herbig-Haro jets 
for example. The maximum asymptotic velocity of a cold super- 
Alfvenic outflow anchored at ro is m Pj00 = Vk,o V2/1 - 3, where 



Table 1. Viscous, accretion, mechanical, kinetic powers and 
MHD Poynting flux. 



Power 



Value 



Viscous Power , P V1SC 0.000936 

Accretion Power, P acc 0.00362 

Mechanical Power, P mec hjet -0.000183 

Kinetic Power, Pkjmaicja 0.000197 

Poynting Flux, P M HD, J(: t 0.000563 



A =i (rA/Vo) 2 is the magnetic lever arm parameter. It is possible 
to estimate A using the ratio: 



F Poynting J s E X B ■ dS 



^kinetic pu 2 u p /2 • dS 



2(A-l) 



(16) 



In our case we derive a A of 2.4, which would provide an asymp- 
totic velocity of 1.4 Vk,o- This is larger than the value of 1.2 
found, which hints to the fact that the magnetic structure retained 
a fraction of its energy. By computing the magnetic contribution 
of the total jet power further out, at z = 100, we found that it still 
represents 33% (Fig.fTOb. This is again in contrast to self-similar 
models where all the MHD power is converted into jet kinetic 
power. This aspect will be discussed in a companion paper. 
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Fig. 3. Time evolution of the ejection rate (top), calculated by 
integrating the mass flux over the super-fast region only, half 
disk accretion rate (middle) measured at the inner radius and the 
ejection to accretion rate ratio (down). After an initial transient 
phase that lasted up to 200 tko, the global system reached a quasi 
steady-state. 



3.2. The SAD structure 

Despite the presence of an outflow (be it a super-FM self- 
confined jet or only super-A flow), the disk structure strongly 
resembles that of a standard accretion disk. 

Most of the released power is radiated away: this implies that 
the main dominant torque is the viscous one. Accretion proceeds 
therefore throughout the disk with a Reynolds number of order 
unity. The value of the accretion Mach number m s = -u r /c s at 



Fig. 4. Poloidal velocity normalized to the Keplerian velocity at 
the footpoint along a magnetic field surface. The time is 953.3 
tko- The radius of the footpoint of the magnetic surface is 2.4. 




Radius, R 



Fig. 5. Ratio of midplane radial accretion sonic Mach number 
m s at time = 953 tko to its theoretical value, a v c s /VK- Accretion 
clearly proceeds at a rate controlled by the anomalous viscosity, 
the jet torque being negligible. 



Z = is a good test as its fiducial value in a S AD should be of the 
order of m± = a v h/r dRozvczka et al.| [l994). Figure [5] shows the 
radial profile of m s /(a v h/r) at a time = 953 tko- The theoretical 
approximation is very close to the simulated one, up to a factor of 
about two. This is very reasonable given the fact that the actual 
expression of the torque involves radial derivatives. The action 
of the jet is therefore totally negligible on midplane accretion. 

Now, let A be the ratio of the magnetic (jet) to the viscous 
torque averaged over the disk thickness. Since most of the avail- 
able power is stored into rotation in a thin accretion disk, the 
fact that the jet carries a tiny fraction of the accretion power is 
directly related to a negligible torque on the bulk of the disk 
mass. Analytically, this may be written 



jet 



A -K 



A 



(17) 



At each radius within zone I, we vertically integrate the torques 
and obtain thereby a distribution A(r). It is relatively smooth, 
with small deviations from an average value of approximately 
0.15. This is consistent with the ratio PMHD.jet/^acc of 0. 1 55 com- 
puted from Table [TJ 

Although the action of the jet is negligible on the equatorial 
accretion motion, it has some impact on the disk vertical struc- 
ture. This is illustrated for instance in Fig. [6] where the density 
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Fig. 6. Vertical density profiles at r — 2.4, 5.5, and r — 14.1, nor- 
malized to their midplane value. At the inner radius, the initial 
steep gradient has been flattened out. This is a signature of mass 
loss from the disk. The time is at 953 r K o. 

profile is shown at different radii. The profiles have been nor- 
malized to the midplane density and plotted against z/h(r) where 
h{r) is the local thermal heightscale. Clearly, the profile located 
at r — 2.4, that is within the ejecting region, becomes flatter than 
that from r = 14.1 (outside the ejecting zone), has a result of 
the mass loss. Notice also the dramatic decrease in density of 
about four decades at the disk surface. We shall come back to 
this feature later on. 
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Fig. 7. Snapshot at t — 953tko of the inner disk regions. The 
color background is the log of mass density with the Fast (dot- 
ted line) and Alfven (dashed line) surfaces overplotted. The last 
super-FM field line is anchored at about r = 5, whereas the last 
super- A field line is at r — 13. The jet exhibits the characteristic 
butterfly shape in the electric current lines (shown in dark blue). 



3.3. Electric currents 

Understanding the behavior of electric currents is the key point 
in accretion-ejection theory. Figure|7]shows a zoom of the eject- 
ing regions I and II at t — 953tko- The Alfven (dashed) and Fast 
(dotted line) critical surfaces are shown along with the poloidal 
electric current lines j p (blue). Globally, some current enters the 
disk at its inner edge (J z < 0) and flows outwardly within the 
disk. In zones I and II where ejection takes place, current lines 
are closed within the jet. The crossing of this poloidal current 
through poloidal field lines results in jet confinement and accel- 
eration. In zone III where there is no jet, thus much less plasma, 
there is almost no current flowing there. 

Let us have a look at the vertical profiles of the magnetic field 
components and electric current density at three radii located 
within the three previously defined zones (Fig. [8]). In all three 
zones, the disk surface can be easily detected as it is the locus 
where (i) the radial component undergoes a huge increase and 
(ii) the toroidal field abruptly changes its behaviour. 

Throughout all the disk, the dominant toroidal current den- 
sity is located at the disk surface not at the disk midplane as 
assumed in the infinitely thin disk approximation. This is easy to 
understand from Ohm's law in resistive MHD. Indeed, neglect- 
ing the contribution of the vertical velocity one gets 



Unless the vertical profile of u r strictly follows that of the mag- 
netic diffusivity, the vertical decrease of v m leads unavoidably 
to an increase of J^. Remarkably, such a profi le of J<t, has been 
alread y discussed in self-similar solutions of ICasse & Ferreiral 
(l2000al) (see their Fig. 7). As a consequence, field lines remain 



straight within the bulk of the disk and bend only at its surface. In 
all zones, such a bending is large enough to satisfy the Blandford 
& Payne energetic criterion for cold wind launching. As a matter 
of fact, the bending of the field lines (B*/B z ) gets larger as one 
goes from zone I to zone III. 

The vertical profile of the toroidal field is controlled by 
the radial current density J r . The profile J r (z) is very interesting: 
contrary to self-similar solutions, most of the poloidal current is 
flowing at the disk surface and not at the disk midplane. Because 
of the small value of B z , the unipolar induction effect is small 
and so is the induced radial current. However, that current be- 
comes much larger towards the disk surface, mostly because of 
the vertical decrease of the resistivity. This results in a magnetic 
shear \Bt/B z \ that goes from 2.25 (zone I) to 6.7 (zone II) or 
even more (> 16 in zone III). A magnetic shear of ~ 2 is a typi- 
cal value already met in previous simulations and in self-similar 
models. It results from the interplay between the disk and the jet 
and leads to a steady-state. Self-similar models have shown that 
a larg er value will result in unsteady disk a nd wind configura- 
tions (iWardle & Koni gill 1 993b iFerreirall 1 997h . Remarkably, zone 
III exhibits an even larger value that still increases with height. 
There is no jet in this zone but only a torsional Alfven wave lead- 
ing to a (negligible) magnetic braking of the underlying disk. In 
this zone, the magnetic field is so small that the shear can be very 
large with no actual damage on the disk structure. 

In zones I and II, where ejection takes place, the radial cur- 
rent density decreases vertically and becomes eventually nega- 
tive. In fact a close examination of Fig. Q reveals the follow- 
ing pattern. The poloidal current enters the disk at its surface 
(J* < 0) between r = 1 and r - 2.5 and then flows inside the 
disk with J, > 0. From r = 2.5 to about r - 13 it exits the 
disk (J z > 0). This is the usual behavior expected in the jet ac- 
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celerating region, with a typical butterfl y shape for the electric 
poloidal current density (see Fig. 13 in Ferreira Hl997l) . This is 
not the case anymore in zone III where no jet is launched. There, 
most of the radial current remains confined within the disk, with 
J* vanishingly small. This is clearly seen in Fig. [7] current lines 
remain within the disk giving rise to a current / = 2nrB,p roughly 
constant with the radius (thus oc r _1 ). Such a radial profile of 
the toroidal field will be discussed in a companion paper. 



3.4. The self-confined jet 

Here we focus on the dynamics of the super-Fast Magnetosonic 
outflow referred to as the jet. Such a structure established on 
a dynamical time scale, namely the local Keplerian time, up to 
r — 5. This is also approximately the time scale for FM waves 
propagating upstream from the FM surface and reaching the disk 
surface. As a consequence, the Keplerian time is also the good 
time scale for establishing a steady-state. Indeed, we observe 
that after roughly 30 Keplerian orbits there is no relevant mod- 
ification of the inner jet structure, corresponding to a few times 
the orbital period at r — 5. 

In order to assess whether or not our magnetized adiabatic 
outflow reached a steady-state, the best way is to compute the 
five following quantities, namely the mass flux to magnetic flux 
ratio: 



k = p 



the specific angular momentum: 
l = ru$- — , 

where = Qr, the magnetic surface rotation rate: 
kB* 

a, = n — - , 

pr 

the entropy: 

and the specific energy or Bernoulli invariant: 



m 2 y P Q*rfi^ 

E= t +$g + ^- t- 2 

2 y — 1 p k 



(19) 



(20) 



(21) 



(22) 



(23) 



According to steady state jet theory these quantities should 
be invariants, namely constant both in time and along each mag- 
netic surface. They are shown along several field lines in Fig. [9] 
as a function of the altitude z at a time 953 tko- All quantities 
first undergo some variation first in the resistive disk region until 
they become constant. The sudden change in their profile (see for 
instance the rise in E and S) occurs at the transition from the re- 
sistive disk to the ideal MHD flow. Further out, it can be verified 
that all quantities are indeed invariants, proving our statement 
that a steady-state outflow has been achieved. 

The specific energy or Bernoulli invariant is separated into 
its kinetic, enthalpy, magnetic and gravitational components (Eq. 
l23l . These are plotted along a single magnetic surface anchored 
at r — 2.4 in Fig. [10] The vertical line at z — 0.9 shows the 
Slow-Magnetosonic point whereas the second line at z - 2 is 
the Alfven point. The vertical line at z - 0.7 shows the alti- 
tude where the resistivity has been set to zero, marking thereby 
the transition between the underlying resistive layers and the 
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Fig. 10. Components of the specific energy E along a magnetic 
surface anchored at r — 2.4: magnetic, kinetic (poloidal and 
toroidal), gravitational and enthalpy. Vertical thick lines indicate 
the heights at which the slow magnetosonic (SM), Alfven (A) 
and fast magnetosonic points (FM) are reached. Notice how neg- 
ligible is the enthalpy compared with other components. 



ideal MHD flow above. Several important aspects can be drawn 
from this plot. First, the enthalpy (solid line) is negligible: 
we are therefore contem plating a "cold" outflow as defined by 
iBlandford & Pavnd d!982l) . The large bending of the magnetic 
field at the disk surface is thereby enough to drive a magneto- 
centrifugal outflow. It can be moreover seen that the dominant 
contribution in E is indeed the magnetic one (+ symbol), as it 
should be in such a case. This magnetic energy is then converted 
into poloidal kinetic energy, but still retains a sizable fraction of 
its initial value at z = 100. This limited efficiency of the energy 
transfer will be discussed elsewhere. 

It is noteworthy that the specific energy E is an invariant only 
after roughly z ~ 1.2: in all the trans-SM zone, it still increases. 
This is not due to the enthalpy as it remains always negligible 
despite the huge increase in entropy at the disk-jet interface (see 
Fig. [9]). This increase in E can be traced back to the increase in 
the magnetic component at that same location (see the increase 
in / in Fig. [9). This is actually due to a decrease of the mass 
to magnetic flux ratio k in the ideal MHD zone. How can this be 
understood in a well tested code where the conservation of quan- 
tities such us mass and total energy is assured up to numerical 
accuracy? 

A numerical algorithm such as the one that we employed 
for our experiments adds to the "ideal" flux a diffusive part 
roughly proportional to the local gradient of the corresponding 
conserved variable, as in the case of the HLL Riemann solver 
used in our simulations. For instance, it can be shown that, in 
a stationary situation, B p ■ Vk is proportional to the divergence 
of a numerical diffusive flux F Pi ns that follows the gradient of 
the density. As a consequence, our estimator of A: is a constant 
only whenever numerical diffusion is really negligible. This is 
clearly not verified at the disk surface where the steepest gradi- 
ent is present. But as we move upwards, the numerical contribu- 
tion vanishes and our estimator converges (decreases) towards 
the real value of k. This points out however to a possible numer- 
ical bias in our MHD simulation. 

Another numerical bias can be related to irreversible numer- 
ical heating, clearly visible in the entropy profiles shown in Fig. 
[9] To test this suspicion, let us have a look at the forces that ac- 



Gareth C. Murphy et al.: Ejection from weakly magnetized disk 



B fields, r= 2.4 



J r , h, r= 2.4 




-0.0010 
-0.0015 



B fields, r= 7.9 




B fields, r=3 1.3 




-2.0'10 




J r , L, r= 7.9 




J r j J*? r— 3 1.3 




Fig. 8. Vertical profiles of the components of the magnetic field (left) and electric current density (right), in zones I (r = 2.4), II 
(r = 7.9) and III (r = 31.3), plotted at a time 953t K o- See text for more details. 



tually drive the poloidal outflow. It is convenient to compute the F m is the parallel component of the Lorentz force 
projection of all forces along a given magnetic surface (Fig.lTTb. 
F v is the parallel component of the kinetic pressure gradient 
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Fig. 9. The five quantities which are considered 
invariant under ideal MHD: the mass flux to 
magnetic flux ratio k; the specific angular mo- 
mentum /; the magnetic surface rotation rate 
fl„; the entropy S and the specific energy or 
Bernoulli invariant E. They are shown in code 
units along field lines anchored (from left to 
right in the plot of 5 ): r = 1 .6, 2, 3, 4. The mea- 
surements are made at t = 953tko- 



where / = 2nrB^ is the electric cu rrent flowing inside that same 
magnetic surface dFerreiralll997l) . F e g is the effective gravita- 
tional + centrifugal force 

B r Q. z r - B ■ Vd> G 

Feff = p j— f , (26) 

\ b p\ 

and, finally, F v is the parallel component of the divergence of the 
shear viscous stress tensor 

g p ■ (V ■ D 

F v = r— : . (27) 

\ B p\ 

Notice that the centrifugal term contained in F e ff arises thanks to 
the azimuthal magnetic force = (J x B) ■ e$ = -B p F m / 'B$. 
With a magnetic shear \B^/B V \ of order unity at the base of the 
jet (Fig. [8}, we have comparable forces F$ ~ F m > 0. 

What is known from analytical studies is that it is mainly 
the vertical component of the plasma pressure gradient that 
lifts the disk material upwards in the r esistive MHD layers^ 
(Ferr eira & Pelletied[l995l lFerreiralll997h . But this effect works 
only in a small vertical extent around the disk surface. This is the 
region where both the radial and vertical velocity components of 
the plasma switch from negative (accretion) to positive (ejec- 
tion). Once this outward movement has been initiated within the 
resistive layers, magnetic and centrifugal forces become both 
dominant and the usual understanding in ideal MHD then ap- 
plies. The critical issue of mass loading, namely the amount of 
mass that is actually ejected from the disk (measured by k), is 
therefore directly related to the delicate interplay of forces in 
this layer. 

This picture is globally confirmed by our numerical simula- 
tion (Fig. [TT1 >. Inside the disk the viscous stress transports mo- 
mentum upwards but it reduces to zero at the disk surface (see 
Appendix [A] for the definition). The projection of the Lorentz 
force is initially negative, showing that the magnetic force within 

3 In the zone where J* < 0, namely where the current enters the disk 
surface (from r = 1 to r - 2.5), the toroidal magnetic pressure provides 
also an upward push (-fl.Sj > 0). But in the zone where Jt > (from 
r — 2.5 to r = 5), the magnetic contribution is only a vertical pinch. 
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Fig. 11. Forces projected along a poloidal magnetic surface an- 
chored at r — 2.4, plotted against the altitude above the disk 
midplane at 953tkq: ^ p is the kinetic pressure gradient, F m the 
Lorentz force, F e ff is the net gravitational+centrifugal forces and 
F v is the viscous stress. The sum of all is also plotted. The ver- 
tical lines indicate the heights where the flow becomes respec- 
tively in ideal MHD (disk surface S d, super-SM {S smX super- A 
(Sa) and super-FM (5fm)- 



the disk and up to its surface is hindering ejection, not helping 
it. The same holds for the effective force, where gravity over- 
comes the centrifugal term even up to after the Alfven point. It 
is indeed the plasma pressure term that makes the difference by 
providing a super-SM ejection. Before the Alfven point however, 
it becomes negligible and the dominant force is the magnetic one 

This enhanced pressure gradient is likely to be related to the 
numerical heating visible in Fig. [9] On the other hand, when 
looking at the Bernoulli invariant, the enthalpy remains always 
negligible: we argue that, due to the enhanced mass flux related 
to the diffusive effects discussed before, the numerical heating 
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Fig. 13. Two simulations done at different resolutions: two times (left panel) and eight times (right panel) our reference simulation. 
The colormap is the log of gas density with overplotted magnetic field lines at footpoints r = 2,3,4,5,6. The overlaid critical 
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Fig. 12. Same as Fig.QT|but with a resolution eight times higher. 
Although the time is now only 15.9tko, a steady-state has been 
already achieved at this radius. 

per unit volume does not correspond to a significant temperature 
and enthalpy increase. 

3.5. Mass loading in numerical simulations 

The side effect of using finite difference methods to solve fluid 
equations is that it introduces numerical biases that play the 
role of a magnetic diffusivity, heat conductivity and viscosity. 
Although such a numerical diffusion is limited so far as is possi- 
ble, it is unavoidable and plays a role wherever a steep gradient 
sets in in any quantity. This is clearly the case for the density 
profile (Pig. O at the resistive-ideal MHD zone where the criti- 
cal is sues of ma s s load ing and initial jet acceleration take place. 
See IZanni et ail (120071) for a discussion of this point. Thus, if 
present, such a numerical diffusion is actually an extra force term 



that will appear in particular in the vertical equation and should 
be present in any numerical simulation published so far. 

The previous clues that some numerical diffusion of mass 
is taking place (the bump in the k "invariant" seen in Fig. [9]) 
can be tested by repeating the simulation at higher resolution. 
We therefore performed two more simulations, one with a res- 
olution twice and another with a resolution 8 times higher (Fig. 
[TJt . The physical parameters, boundary and initial conditions re- 
mained unchanged, but the physical extent of the simulation was 
reduced. Additionally the simulations were only carried out for 
~ 17 inner disk orbital periods. 

We obtained the following results : (i) whereas a super-FM 
jet is still launched in a steady-state, (ii) the radial extent of the 
ejecting zone is narrower up to r ^ 2 only. We will come to back 
to this later. Figure [14] shows the various quantities that should 
remain invariant, to be compared with Fig.|9]with the lowest res- 
olution. The anchoring radii are the same as in Fig. [9] Clearly, the 
invariants are flatter as numerical diffusion is reduced. Also, the 
bump in k has now almost vanished and the transition from resis- 
tive to ideal MHD is much better caught. Moreover, the entropy 
profiles clearly shows that the numerical heating is strongly re- 
duced. 

Figure Q~2] shows the parallel forces along a given field line 
anchored at the same radius as in the lowest resolution case. 
The general trend remains the same although the effect of the 
thermal push is now dramatically reduced. It is still the domi- 
nant force allowing a trans-SM flow but its importance decreases 
more rapidly. Of course, only magnetic forces provide a super- 
Alfvenic flow. The reduction of the ejection efficiency with in- 
creasing resolution confirms our suspicion: numerical diffusion 
is indeed at work at the disk surface in the inner regions of the 
grid. This effect naturally explains the mass loading, initial push 
and thereby increase in the specific energy E. 
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Fig. 14. The same as Fig. [9] but for a resolu- 
tion eight times higher. The invariants are flat- 
ter and the bump in the mass loading is reduced. 
A steady-state super-FM jet is still present but 
from a smaller radial extent. The measurements 
are at time 17.5tko- 



3.6. The role of the disk magnetization 

So far we understood the initial mass loading and driving mech- 
anisms of the jet. However, what determines its radial extent 
has remained unexplored. In previous simulations of that kind 
the extent of the ejecting zon e was increasing in time with 
a Keplerian scal ing r oc t 2/3 dCasse & Keppensll200l [20041 
Zannietal . 2007) . It is the first time where a jet is launched from 
a finite region that remains constant over time. 

In steady-state jet theory, the Bernoulli invariant must be 
positive at all magnetic surfaces. Neglecting enthalpy (as Fig. 
[TO] suggests), Eq. ( 1231 provides 



o n cr — 1 

F ~ Ci 2 r 2 

c _ "*,()'() 2 

where 



(28) 



cr - — 



pu 2 u v 



(29) 



z=h 



is the ratio of the MHD Poynting flux to the kinetic energy flux 
measured at the disk surface. This quantity is sometimes referred 
to as the (initial) jet magnetization. A cold jet requires therefore 
cr + larger than unity. We plotted in Fig.[T5l(top) the jet magneti- 
zation as function of the disk launching radius for our reference 
simulation. Beyond a radius of about 5, this quantity becomes in- 
deed smaller than unity, corresponding nicely to the end of zone 
I (super-FM jet). A super-A outflow is nevertheless launched at 
larger radii, but this is a matter dominated flow that never reaches 
a steady-state. The overall picture is therefore consistent. But 
what determines the radial distribution cr + (r)? 

Another way to write the initial jet magnetization is 



<r + = 2f— *- 

f u ' 



(30) 



where = B 2 JP + is measured at the disk surface. We plotted 
in Fig. [15] (down) the disk magnetization at the upper surface 
layers. It can be seen that cr + (r) follows approximately the same 
trend as ju + (r), namely a radial decrease. 



Fig. 15. Upper panel: Initial jet magnetization cr + (Poynting to 
kinetic flux ratio) measured at the disk surface in the inner 
zones of the accretion disk. Lower panel: Disk magnetization 
= B 2 /P + ) measured at the disk surface. These curves were 
obtained with our reference simulation at the final stage. 



In fact, analytical calculations done within the self-similar 
framework already pointed out the importance of the disk mag- 
netization y u for launching super-FM jets. It was sh own that 
isothe rmal dFerreira & Pelletierl 119951: iFerreiral 1 19971) or adia- 
batic (ICasse & Ferr eira 2000a) magnetic surfaces require a field 
close to equipartition, namely fi smaller but around unity. Our 
own results suggest that it is the disk magnetization that actu- 
ally defines the ejecting zones. Beyond r — 5, the magnetic field 
would be too small to allow a proper jet to be launched. 

Let us make a very crude approach by assuming a B z com- 
ponent almost constant in the vertical direction and an isother- 
mal hydrostatic density profile. In that case, one would have 
/i(z) = yuexp(z 2 /2h 2 ) where fi is the disk magnetization at the 
disk midplane. It seems therefore dubious that yu + , reaching a 
value 10 to 100 times /j at a few scale height, could ever reach 
a value of order unity if /j is too small. So, even in the presence 
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Fig. 16. The colourmap shows the logarithm of gas density 
in the weak magnetic field simulation at time t=321.5rKo. The 
dashed lines enclose material moving at velocities faster than the 
local Alfven speed. There is no superfast outflow present and the 
superAlfvenic material is extremely fragmented. 



of a numerical diffusion, no jets should be produced if the disk 
magnetization is too low. 

In order to test this conjecture, we performed another nu- 
merical simulation with fi decreased by one order of magnitude 
(namely starting at 2 x 10~ 4 at the disk inner radius). The other 
physical parameters, boundary and initial conditions, as well as 
the numerical resolution were otherwise identical to the refer- 
ence simulation. We found that, in this case the super-Alfvenic 
material was extremely sporadic and fragmented in the domain, 
and no super-fast-magnetosonic jet was observed. See Fig. [16] 

In our view, this clearly confirms that the disk magnetization 
must be high enough in order to launch self-confined (super-FM) 
jets. This result goes in the same direction as those obtained with 
self-similar solutions. However, the latter claimed that only /j. 
smaller but close to unity (namely a field close to equipartition) 
allows the launching of magnetized jets. The physical argument 
is the following. For a jet to be launched, the lifted mass must 
cross the SM point around the disk surface. In a cold environ- 
ment the only force able to do this is the magnetic one. It turns 
out however, that it is much easier to do it if the accretion ve- 
locity is already not too far from the sound speed. This is the 
reason why isothermal or adiabatic jets require fields close to 
equipartition (Ferreira & Casse, 2009, submitted). 

Apparently, this is in contradiction with our own result since 
we do obtain jets with /i of the order of a few 10~ 3 . The rea- 
son for this discrepancy lies in the fact that the analytical mod- 
els were obtained under the assumption of either isothermal or 
adiabatic magnetic surfaces. Here, as we showed, there is a nu- 
merical diffusion that allowed mass to leak from the disk to the 
open, rotating field lines. This extra eff ect has been mimicked 
for instance in ICasse & Ferreira! (l2000bl) with the presence of a 
heating term at the disk surface. New solutions, called "warm" 
in contrast to the previous "cold" ones, were found with an en- 



hanced mass flux. But these authors did not recognize that the 
required disk magnetization [i was indeed smaller than for cold 
jets. We report here that it is indeed the case, with some of the 
"warm" self-similar solutions found with fi = 0.08. 



4. Concluding remarks 

In this paper, we performed four 2.5D numerical MHD simu- 
lations of a resistive viscous accretion disk threaded by a weak 
magnetic field. The initial magnetic field distribution was chosen 
so that the disk magnetization /i = decreases radially from 
the central object. Our reference simulation, done with a maxi- 
mum value of fi — 2 x 10~ 3 , has run for more than 950 Keplerian 
orbits at the inner radius and is therefore the longest to date. 

It is shown that the disk structure resembles that of a stan- 
dard Shakura & Sunyaev disk with accretion controlled by the 
turbulent (alpha) viscous torque only. However, a super fast mag- 
netosonic, self-confined jet is observed to be launched from the 
inner disk regions. It is first time that (i) steady-state super-FM 
jets are produced from a weakly magnetized disk and (ii) from 
a finite disk region that remained constant over time. The power 
carried away by these jets is tiny and directly related to the neg- 
ligible torque on the disk. The dynamics of the jet and its prop- 
agation into the medium will be studied in a forthcoming paper. 
Here, we focused on the jet acceleration region where the flow 
crossed the three MHD critical surfaces (Slow Magnetosonic, 
Alfven and Fast Magnetosonic). 

The critical issues of mass loading and initial jet accelera- 
tion (the crossing of the SM surface) are shown to be strongly 
affected by the unavoidable steep decrease of the density pro- 
file at the disk surface. Such an effect has been underestimated 
in previous simulations. It is the quality of the grid resolution at 
the disk surface that ultimately determines the amount of ejected 
mass. One way to solve this problem is to use either an enhanced 
resolution at the disk surface, a less diffusive algorithm, a higher 
order method or an adaptive grid which refines on the density 
gradient. 

We argue however that this feature might mimic some 
additional h eat input at the disk s urface, as explo r ed fo r 
instance by iQgilyie & Livid d 19981) . lOgilvie & Liyiol d200lh . 
ICasse & Ferreira! (I2000b). This aspect is extremely promising as 
most astrophysical accretion disks do probably have superheated 
layers due to irradiation by the central object (young stars, cata- 
clysmic variables) and/or some X-ray source (e.g. around black 
holes). As a consequence, "cold" (e.g. isothermal or adiabatic) 
ejection is probably never achieved in Nature. 

This allows also to relax the constraint of equipartition fields 
needed for driving jets as our jets were obtained from a very low 
magnetized disk (but not too low). This opens a new fascinat- 
ing topic: the magnetic history of any given object. One might 
indeed consider accretion disks displaying a whole continuum 
in ejection efficiency, from jets carrying a sizable fraction (if 
not most) of the released accretion power to jets that are a mere 
epiphenomenon of accretion. For any given object, the key pa- 
rameter would be the disk magnetization. This clearly deserves 
further investigation. 
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where Pdo = e Pd()V^ (y The components of the poloidal speed 
are: 
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The toroidal speed is: 
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For viscosity and r esistivity, the expression used by 
IZanni & Ferreira! ((2009) is employed: 
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where the isothermal soundspeed calculated on the disk mid- 
plane c s (r)\ z= o can change in time. 

For the atmosphere above disk all velocities are set to zero, 
u r = u z = u,p = 0, and a hydrostatic, spherically symmetric 
atmosphere is prescribed: 
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A density contrast p a o/Pdo = 10 4 has been assumed in all 
the simulations. _ 

The components of the viscous stress tensor T used in 
PLUTO are: 



du r ( 2 , 

T rr = 2t 1v-7^- + U - ~77v I V ■ H 



(A.9) 



Appendix A: Additional numerical conditions 

As an initial condition for the simulation the perturbative solu- 
tion of the steady-state MHD equations is taken. The disk pres- 
sure and density are computed by solving the hydrostatic verti- 
cal equilibrium, the toroidal speed is determined by the radial 
equilibrium, whereas the radial velocity is given by the angular 
momentum conservation equation. We assumed a thermal disk 
heightscale h = er, where the aspect ratio e = c s /Vk = 0.1 
is given by the ratio between the isothermal soundspeed c s = 
■JP/p and the Keplerian speed Vk 



yGM/r calculated at the 
disk midplane. The disk density is therefore given by: 



Pd = PdO • 
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and the thermal pressure by: 



P* = PJ^) 
\Pd0j 



(A.2) 



where the bulk viscosity £ is set to zero. 



